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Abstract. Self-induced flavor conversion of supernova (SN) neutrinos is a generic feature 
of neutrino-neutrino dispersion. The corresponding run-away modes in flavor space can 
spontaneously break the original symmetries of the neutrino flux and in particular can spon¬ 
taneously produce small-scale features as shown in recent schematic studies. However, the 
unavoidable “multi-angle matter effect” shifts these small-scale instabilities into regions of 
matter and neutrino density which are not encountered on the way out from a SN. The tra¬ 
ditional modes which are uniform on the largest scales are most prone for instabilities and 
thus provide the most sensitive test for the appearance of self-induced flavor conversion. As 
a by-product we clarify the relation between the time evolution of an expanding neutrino 
gas and the radial evolution of a stationary SN neutrino flux. Our results depend on several 
simplifying assumptions, notably stationarity of the solution, the absence of a “backward” 
neutrino flux caused by residual scattering, and global spherical symmetry of emission. 
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1 Introduction 

Core-collapse supernovae (SNe) or neutron-star mergers are powerful neutrino sources and 
probably the only astrophysical phenomena where these elusive particles are dynamically 
important and crucial for nucleosynthesis [1, 2]. The low energies of some tens of MeV imply 
that /3 reactions of the type v e + n «->• p + e~ and u e + p ^ n + e + are the dominant charged- 
current processes. Heavy-lepton neutrinos 7^, v T1 and P T , in this context often collectively 
referred to as v x , interact only by neutral-current processes. Therefore, neutrino energy 
transfer and the emitted fluxes depend on flavor and one may think that flavor oscillations 
are an important ingredient. However, the large matter effect in this dense environment 
implies that eigenstates of propagation and those of interaction very nearly coincide [3, 4], 
In spite of large mixing angles, flavor oscillations are irrelevant except for MSW conversion 
when neutrinos pass the resonant density as they stream away [5, 6]. Therefore, the neutrino 
signal from the next galactic SN may carry a detectable imprint of the yet unknown neutrino 
mass hierarchy [7-11]. 

This picture can fundamentally change when the refractive effect of neutrinos on each 
other is included [12-14]. The mean field representing background neutrinos can have fla¬ 
vor off-diagonal elements (“off-diagonal refractive index”) due to flavor coherence caused by 
oscillations and can lead to strong flavor conversion effects [15-54]. (We ignore additional 
effects that would arise from non-standard neutrino interactions [55], spin-flip effects caused 
by neutrino magnetic dipole moments [56-58], by refraction in inhomogeneous or anisotropic 
media [59-62], or the role of neutrino-antineutrino pair correlations [63-66].) Self-induced 
flavor conversion preserves the global flavor content of the ensemble, but re-shuffles it among 
momentum modes or between neutrinos and antineutrinos. The simplest example would be 
a gas of v e and v e converting to and 7^, leaving the overall flavor content unchanged. 
The interacting modes of the neutrino field can be seen as a collection of coupled oscillators 
in flavor space. The eigenmodes of this interacting system include collective harmonic os¬ 
cillations, but can also include run-away solutions (instabilities) which lead to self-induced 
flavor conversion [32]. Under which physical conditions will instabilities occur, how can we 
visualize them, and how will they affect the flavor composition of neutrinos propagating in 
the early universe or stream away from a SN core? 

To study these questions, many simplifications were used and especially symmetry as¬ 
sumptions were made to reduce the dimensionality of the problem. However, symmetry 
assumptions suppress those unstable solutions which break the assumed symmetry. There¬ 
fore, when instabilities are the defining feature of the dynamics, symmetry assumptions 
about the solutions can lead to misleading conclusions because, even if the system was set 
up in a symmetric state, the interacting ensemble can break this symmetry spontaneously. 
This behavior is analogous to the hydrodynamical aspects of SN physics which cannot show 
convective overturn if the simulation is spherically symmetric, yet such 3D effects are now 
understood to be crucial for SN physics [67-82] . 

Our present concern is the question of “spatial spontaneous symmetry breaking” in self- 
induced flavor conversion. In previous studies, the flavor content of neutrinos streaming away 
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from a SN core was taken to remain uniform in the transverse directions. However, recent 
studies of simplified systems suggest that this symmetry can be spontaneously broken [51— 
54]. To avoid the complication of global spherical coordinates, it is sufficient to model the 
neutrino stream at some distance with plane-parallel geometry, i.e., we can use wave vectors 
k in the transverse plane to describe small-scale spatial variations. In this terminology, 
traditional studies only considered k = 0 (global spherical symmetry). In analogy to this 
k = 0 case it was found that for any k there is some range of effective neutrino densities 
where unstable solutions exist [51-54]. We usually express the neutrino density in terms of 
an effective neutrino-neutrino interaction energy ji = \/2G-Fn Ue {r) (R/r) 2 , where n Ve (r) is 
the v e density at distance r and R is some reference radius playing the role of the neutrino 
sphere. The parameter /r varies with r~ 4 because the neutrino density decreases as r~ 2 with 
distance. In this terminology, for any k there is a range /i m j n < fi < /r max where the system is 
unstable. For larger k (smaller spatial scales), the instability range shifts to larger /r, i.e., to 
regions closer to the SN core. This finding suggests that the neutrino stream is never stable 
because at any neutrino density there is some range of unstable k modes. 

However, such conclusions may be premature as one also needs to include the refractive 
effect of matter which also tends to shift the instability to larger /r, a phenomenon termed 
“multi-angle matter suppression” of the instability [26, 37]. We usually express the multi¬ 
angle matter effect in terms of the parameter A = \/2G-pn e {r ) ( R/r ) 2 , where n e (r ) is the 
electron density at distance r. One should study the instability region in the two-parameter 
space of effective matter and neutrino densities, A and pi, which we call the “footprint of the 
instability.” In figure 1 we show as an example the footprint of the MAA instability for a 
schematic SN model (MAA stands for “multi azimuth angle,” i.e., one type of instability). 
We always define the instability region by the requirement that the growth rate k > IIP 2 uo, 
where ujq is a typical vacuum oscillation frequency. 

Our schematic SN model consists of neutrinos and antineutrinos with a single energy 
corresponding to a vacuum oscillation frequency cuo = Am 2 /2E = 1 km -1 . We assume they 
are emitted isotropically at the neutrino-sphere radius R = 30 km. We consider two-flavor 
oscillations and, after subtracting the v x flux, we are left with a v e and v e flux such that there 
emerge twice as many v e than v e . Finally, we assume that the effective neutrino-neutrino 
interaction energy, n = \/2GFn Ue (R/r) 2 , at the neutrino sphere r = R is /Uo = 10 5 km -1 . 
Under these circumstances, run-away solutions for k = 0 exist within the footprint area 
shown in figure 1 as a blue-shaded region. We also show, as a solid red line, a possible density 
profile, corresponding to electron density as function of radius. The sudden density drop at 
r ~ 200 km is the shock wave. In this example, the density profile does not intersect with 
the instability footprint for radii below the shock wave. On the other hand, for larger radii, 
collective flavor conversion would begin. We further show the footprint for inhomogeneities 
with assumed wave-number k = 10 2 and k = 10 3 , measured in units of the vacuum oscillation 
frequency ujq. Notice that we define k in “co-moving” coordinates along the radius, i.e., a 
fixed k represents a fixed angular scale, not a fixed length scale. Notice also that for the 
MAA instability shown here, which is relevant for normal mass ordering, non-zero k values 
lead to two instability regions. This phenomenon does not arise for the bimodal instability. 

The full range of all fc-values inevitably fills the entire region below the k = 0 mode 
(blue shading) in this plot, i.e., the entire gray-shaded region is unstable, whereas the region 
above the blue-shaded part remains unscathed. In other words, essentially the largest-scale 
mode with k = 0 is the “most dangerous” mode. If this mode is stable on the locus of the SN 
density profile in figure 1, the higher-A; modes are stable as well. Of course, if any instability 
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Figure 1. Footprint of the MAA instability region in the parameter space of effective neutrino density 
fj, = \/2G-Fn V ',(R/r ) 2 , where R is the neutrino-sphere radius, and matter density A = V2GFn c (R/r) 2 
for the schematic SN model described in the text. Because p, ex r -4 , the horizontal axis is equivalent 
to the distance from the SN as indicated on the lower horizontal axis. We also show a representative 
schematic SN density profile where the sharp density drop marks the shock wave. We also show 
the instability footprint explicitly for co-moving wave numbers k = 10 2 and k = 10 3 in units of the 
vacuum oscillation frequency. Notice that for the same value of k there are two separate instability 
strips. The collection of all small-scale instabilities fill the gray-shaded region below the traditional 
k = 0 (blue shaded) instability region, whereas they leave the space above untouched. 


is encountered by the physical SN density profile, these instabilities will span a range of scales 
and create complicated flavor conversion patterns. 

The rest of our paper is devoted to substantiating this main point and to explain our 
exact assumptions. We stress that our simplifications may be too restrictive to provide a 
reliable proxy for a realistic SN. In particular, we assume stationary neutrino emission and 
that the solution is stationary as well, i.e., we assume that the evolution can be expressed 
as a function of distance from the surface alone. We also ignore the “halo flux” caused by 
residual scattering which can be a strong effect. Our study would not be applicable at all 
in regions of strong scattering, i.e., below the neutrino sphere. We assume that the original 
neutrino flux is homogeneous and isotropic in the transverse directions, i.e., global spherical 
symmetry of emission at the neutrino sphere. It has not yet been studied if this particular 
assumption has any strong impact on the stability question, i.e., if violations of such an 
ideal initial state substantially change the instability footprint, or if such disturbances would 
simply provide seeds for instabilities to grow. It is impossible to understand and study all 
effects at once, so here we only attempt to get a grasp of the differential impact of including 
spatial inhomogeneities in the form of self-induced small-scale flavor instabilities. All the 
other questions must be left for future studies. 
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2 Equations of motion 


Beginning from the full equation of motion for the neutrino density matrices in flavor space, 
we develop step-by-step the simplified equations used in our linearized stability analysis. In 
particular, we formulate the stationary, spherical SN problem, where the flavor evolution is 
a function of radius, as an equivalent time-dependent 2D problem in the tangential plane. 
A fixed neutrino speed in the tangential plane corresponds to the traditional “single angle” 
treatment, whereas neutrino speeds taking on values between 0 and a maximum, determined 
by the distance from the neutrino sphere, corresponds to the traditional “multi angle” case. 


2.1 Setting up the system 

We describe the neutrino field in the usual way by 3x3 flavor matrices g(t, r, E, v), where the 
diagonal elements are occupation numbers for the different flavors, whereas the off-diagonal 
elements contain correlations among different flavor states of equal momentum. We follow 
the convention where antineutrinos are described by negative energy E and the corresponding 
matrix includes a minus sign, i.e., it is a matrix of negative occupation numbers. 

We always work in the free-streaming limit, ignoring neutrino collisions. In this case, 
neutrino propagation is described by the commutator equation [14, 45] 

i(a t + v-V r )e=[H,g], (2.1) 


where g and H are functions of t, r, E, and v. The Hamiltonian matrix is 




r ( v - v ')2 

+ / dF' --- Qt,r,E'y 


( 2 . 2 ) 


where M 2 , the matrix of neutrino mass-squares, is what causes vacuum oscillations. The 
matrix of charged-lepton densities, N^, provides the usual Wolfenstein matter effect. The 
integration dT' is over the neutrino and antineutrino phase space. Because antineutrinos are 
denoted with negative energies, we have explicitly f dT’ = dE'E' 2 f dv’ / (2-7r) 3 and the 
velocity integration dv' is over the unit sphere. Because the neutrino speed |v| = 1 we were 
able, for later convenience, to write the current-current velocity factor in the unusual form 

(1 —vV) = i(v-v') 2 . 

Studying this 7-dimensional problem requires significant simplifications. For neutrino 
oscillations in the early universe, one will usually assume initial conditions at some time 
t = 0 and then solve these equations as a function of time. To include spatial variations, one 
may Fourier transform these equations in space, replacing the spatial dependence on r by a 
wave-number dependence k, whereas v • V r —> iv ■ k and the r.h.s. becomes a convolution 
of Fourier modes [51]. One can then perform a linearized stability analysis for every mode 
k and identify when modes of different wave number are unstable and lead to self-induced 
flavor conversion [52]. One can also use this representation for numerical studies [51, 53, 54], 

The other relatively simple case is inspired by neutrinos streaming from a supernova 
(SN) core. One assumes that, on the relevant time scales, the source is stationary and that the 
solution is stationary as well, so that dt —> 0. In addition, one assumes that neutrinos stream 
only away from the SN so that it makes sense to ask about the variation of the neutrino 
flavor content as a function of distance, assuming we are provided with boundary conditions 
at some radius R which we may call the neutrino sphere. Actually, this description can be 
a poor proxy for a real SN because the small “backward” flux caused by residual neutrino 
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scattering in the outer SN layers, the “halo flux,” can be surprisingly important for neutrino- 
neutrino refraction because of its broad angular range [40-42] . Here we will ignore this issue 
and use the simple picture of neutrinos streaming only outward. 

We stress that this simplification is the main limitation of our study and its interpre¬ 
tation in the physical SN context. If (self-induced) instabilities exist on small spatial scales, 
they could even exist below the neutrino sphere where the picture of neutrinos streaming only 
in one direction would be very poor. The approach taken here to reduce the 7-dimensional 
problem to a manageable scope may then hide the crucial physics. Therefore, our case study 
leaves open important questions about a real SN. 

2.2 Large-distance approximation 

The main point of our study is to drop the assumption of spatial uniformity, i.e., we include 
variations transverse to the radial direction. However, we are not interested in an exact 
description of large-scale modes. At some distance from the SN, outward-streaming neutrinos 
cannot communicate with others which travel in some completely different direction as long 
as we only include neutrino-neutrino refraction and not, for example, lateral communication 
by hydrodynamical effects. If we are only interested in relatively small transverse scales, 
we may approximate a given spherical shell locally as a plane, allowing us to use Cartesian 
coordinates in the transverse direction rather than a global expansion in spherical harmonics. 

We now denote with z the “radial” direction, and use bold-faced characters to denote 
vectors in the transverse plane, notably a for the coordinate vector in the transverse plane 
and (3 the transverse velocity vector. In the stationary limit, the equation of motion (EoM) 
becomes 

i(v z d z + (3-V a )g=[H,g\, (2.3) 

where v z = \/1 — /3 2 and g and H depend on z, a, E, and (3. If the neutrino-sphere radius 
is R, then at distance r from the SN the maximum neutrino transverse velocity is /3 max ~ 
R/r <C 1. This latter “large distance approximation” is the very justification for using 
Cartesian coordinates in the transverse direction. 

Therefore, it is self-consistent to expand the equations to order /3 2 . (We need to go to 
quadratic order lest the neutrino-neutrino interaction term vanishes.) The /3 expansion is 
not necessary for our stability analysis, but avoiding it does not provide additional precision 
and performing it provides significant conceptual clarity. Numerical precision for a specific 
SN model is not our goal, and in this case we would have to avoid modeling the transverse 
direction as a flat space anyway, especially when considering regions that are not very far 
away from the SN core. 

In equation (2.3) we multiply with l/v z ~ (1 + \fi 2 ) and notice that the gradient term 
remains unchanged if we expand only to order f3 2 so that 

i(d z +(3-V a )g = [H,g], (2.4) 

where the Hamiltonian matrix is the old one times (1 + \(3 2 ) or explicitly 

H = i 1 + ?) ( 2 ^ + / dT ' (/3 ~f )2 ez,*,E',p> ■ (2.5) 

The flux factor under the integral in the second expression is also an expansion to 0(/3 2 ) in 
the form 1 - v z v' z - (3 • (3' = 1 - y/l - - /3' 2 - (3 ■ f3' « \p 2 + ±/3' 2 - (3 ■ (3' = \{(3 - (3') 2 . 

Multiplying this expression with (1 + \t3 2 ) makes no difference because it is already 0(/3 2 ). 
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As a next step, we re-label our variables and denote the radial direction z as time t. 
Moreover, we rescale the transverse velocities as (3 = v/3 max where v is now a 2D vector 
obeying 0 < |v| < 1. Coordinate vectors in the transverse plane are also rescaled as x = 
a/3max 5 he., the new transverse coordinate vector x is “co-moving” in that it denotes a fixed 
angular scale relative to the SN. After these substitutions, the EoMs are 

i{d t + v ■ V x )£> = [H, g] , (2.6) 

with 

H = ^1 + V 2 ^ + \/2GfN^ + V^GF^max J ( ^' - -^~ Qt^,E',v' • (2.7) 

Of course, the neutrino phase-space integration f dr' is understood in the new variables. 

Our stationary 3D problem has now become a time-dependent 2D problem. In the SN 
interpretation, the aspect ratio of the neutrino sphere shrinks with distance and correspond¬ 
ingly, f3 max shrinks. In other words, the physics is analogous to neutrino oscillations in the 
expanding universe. In the SN case, linear transverse scales grow as r/R , where r is the dis¬ 
tance to the SN, i.e., our “Hubble parameter” is R where R is the neutrino-sphere radius 
and the scale factor grows linearly with “time.” The physical neutrino density decreases with 
inverse-distance squared and in addition, the factor /3^ ax accounts for the decreasing value 
of the current-current factor in the neutrino interaction. Therefore, the effective neutrino 
number density decreases with (scale factor)in the familiar way. 

2.3 Single angle vs. multi angle 

In the SN context, one often distinguishes between the single-angle and multi-angle cases, 
referring to the zenith angle of neutrino emission at the SN core. If all neutrinos were 
emitted with a fixed zenith angle, their transverse speeds would be \f3\ = /3 max and in our 
new variables, |v| = 1. In this case (1 + |/?max v2 ) -> (1 + is simply a small and 

negligible numerical correction to the vacuum oscillation frequencies and the matter effect. 
Also, we can revert to the traditional form of the flux factor ^(v — v 7 ) 2 = (1 — v ■ v 7 ). 
Therefore, the SN single-angle case is equivalent, without restrictions, to a 2D neutrino gas 
evolving in time. Therefore, neutrino oscillations in an expanding space (“early universe”) is 
exactly equivalent to the single-angle approximation of neutrinos streaming from a SN core 
with properly scaled effective neutrino and matter densities. 

It has been recognized a long time ago that in the single-angle case, the ordinary matter 
effect has no strong impact on self-induced flavor conversion [19]. As usual, one can go to a 
rotating coordinate system in flavor space. In this new frame, the matrix of vacuum oscillation 
frequencies, M 2 /2 E, has fast-oscillating off-diagonal elements and, in a time-averaged sense, 
it is diagonal in the weak-interaction basis. These fast-oscillating terms are what kick-starts 
the instabilities at the beginning of self-induced flavor conversion but are otherwise irrelevant. 
For a larger matter effect, more e-foldings of exponential growth of the instability are needed 
to “go nonlinear.” In this sense, matter has a similar effect concerning the onset of the 
instability that would be caused by reducing the mixing angle. These effects concern the 
perturbations which cause the onset of instabilities, not the existence and properties of the 
unstable modes themselves. 

We may ignore the small correction to the vacuum oscillation frequency provided by 
the factor (1 + g/j ^, ax v 2 ). We need to keep terms of order /3 2 in the context of the matter 
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and neutrino-neutrino term which in the interesting case are large and, after multiplication 
with /3 2 iax , still larger than the vacuum oscillation term. Therefore, we find 

H = ( + v^GVN^Lx y + V2G F (3 2 max J dV (V ~ v/)2 Q t ^E>y , (2.8) 

where the first term symbolizes the time-averaged vacuum term in the fast co-rotating frame. 
In the single-angle case, where v 2 = 1 for all modes, the remaining matter term can be rotated 
away as well. 

The multi-angle SN case, in this representation, corresponds to a 2D neutrino gas with 
variable propagation speed 0 < |v| < 1 , i.e., the velocity phase space is not just the surface 
of the 2D unit sphere (a circle with |v| = 1), but fills the entire 2D unit sphere (a disk with 
| v| < 1). There is no counterpart to this effect in a “normal” neutrino gas. The early-universe 
analogy does not produce multi-angle effects, but we can include them without much ado by 
allowing the neutrino velocities to fill the 2D unit sphere. 

If neutrinos are emitted “black-body like” from a spherical surface, from a distance this 
neutrino sphere looks like a disk of uniform surface brightness, in analogy to the solar disk in 
the sky. Therefore, this assumption corresponds to the neutrino transverse velocities filling 
the 2D unit sphere uniformly. In earlier papers of our group, we have used the variable u = v 2 
with 0 < u < 1 as a co-moving transverse velocity coordinate, representing the neutrino 
zenith angle of emission. In terms of this variable, the black-body like case corresponds to 
the familiar top-hat u distribution on the interval 0 < u < 1 . 

Our overall set-up was inspired by that of Duan and Shalgar [52], except that they 
consider only one transverse dimension with single-angle emission at the SN. In other words, 
their system is equivalent to two colliding beams evolving in time and allowing spatial vari¬ 
ations. A numerical study of this case, in both single and multi angle, was very recently 
performed by Mirizzi, Mangano and Saviano [53], going beyond the linearized case. Not un¬ 
expectedly, the outcome of the nonlinear evolution is found to be flavor decoherence. In our 
approach, multi-angle effects in this colliding-beam system can be easily included by using 
a ID velocity distribution that fills the entire interval — 1 < v < +1 and not just the two 
values v = ±1. A black-body like zenith-angle distribution corresponds to a uniform velocity 
distribution on this interval. 


2.4 Two-flavor system 

As a further simplification we limit our discussion to a two-flavor system consisting of v e and 
some combination of and v T that we call i / x , following the usual convention in SN physics. 
We are having in mind oscillations driven by the atmospheric neutrino mass difference and 
by the small mixing angle © 13 . The vacuum oscillation frequency is 

10 MeV\ 

E ) 


(2.9) 


LJj = 


Am 2 

2E 


= 0.63 km 


-1 


where we have used Am 2 = 2.5 x 10 ~ 3 eV 2 . Henceforth we will describe the neutrino energy 
spectrum by an 00 spectrum instead, with negative 00 describing antineutrinos. 

The matrix of vacuum oscillation frequencies, in the fast-rotating flavor basis, takes on 
the diagonal form 



( 2 . 10 ) 
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where we have removed the part proportional to the unit matrix which drops out of com¬ 
mutator expressions. We do not include the fast-oscillating off-diagonal elements which is 
irrelevant for the stability analysis. The matter effect appears in a similar form, 


V2G F N e 0l^ 


Av 2 

2 



( 2 . 11 ) 


where again we have removed the piece proportional to the unit matrix. The parameter 
A describing the multi-angle matter effect at distance r from the SN with neutrino-sphere 
radius R and using /3 max = R/r is 

A = y/2Gprieir) 1 ^ = 3.86 x 10 8 km' 1 ^ 12 ^ ~2 > (2.12) 

r z 10 iZ g cm' ,i r z 

where n e (r) is the net density of electrons minus positrons, p(r) the mass density, and Y e (r) 
the electron fraction per baryon, each at radius r. The matter density drops steeply outside 
the neutrino sphere and jumps downward by an order of magnitude at the shock-wave radius. 
Therefore, we need to consider A values perhaps as large as some 10' km -1 all the way to 
vanishingly small values. 

Turning to the neutrino-neutrino term, notice that the g matrices play the role of 
occupation numbers and that the f dT integration includes the entire phase space of occupied 
neutrino and antineutrino modes. Therefore, = J dTg is a flavor matrix of net neutrino 
minus antineutrino number densities, in analogy to the corresponding charged-lepton matrix 
N/ 1 . It is less obvious, however, how to best define an effective neutrino-neutrino interaction 
strength p which plays an analogous role to A. If we were to study a system that initially 
consists of equal number densities of v e and P e , the matrix vanishes, but later develops 
off-diagonal elements. Therefore, we rather use the number density of v e without subtracting 
the antineutrinos and define 


p = y/2Gyn Ue (r) 


R 2 


4.72 x 10 5 km -1 


L Ue 10 MeV 
4xl0 52 erg/s (E Ve ) 


30 km 
R 



(2.13) 


where L Ve is the v e luminosity and (E Ue ) their average energy. More precisely, n Ve is the v e 
density at radius r that we would obtain in the absence of flavor conversions after emission 
at radius R. Previously we have sometimes normalized p to np e instead, or to the difference 
between the P e and P x densities. However, in our schematic studies we assume that initially 
we have only a gas consisting of v e and P e , again obviating the need for these fine distinctions. 
The exact definition of p has no physical impact because it always appears as a product with 
the density matrices. 

In previous papers [32, 45], a further factor 1/2 was included in the definition of the 
multi-angle A and p. We have kept this factor explicitly in equation (2.8) both in the matter 
term and in the flux factor ^(v — v 7 ) 2 to maintain its traditional form. In this way, the 
equations can be directly applied to a traditional “early universe” system. To make contact 
with previous SN discussions, one can always absorb this factor in the definition of A and p. 

As a next step, we project out the trace-free part of the density matrices and normalize 
them to account for the above normalization of the effective neutrino-neutrino interaction 
strength p. 


£?t,x,o;,v 


T r (ft,x,cj,v) , n v e r 

2 ' 2 ' J ^5 X 5 Ct, 5 V * 


(2.14) 
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With these definitions, the two-flavor EoMs finally become 


V * ^x)Gt ? x,a;,v — [Ht } x,o;,V) Gf,x,cj,v] , 


(2.15) 


with the Hamiltonian matrix 


H 


£,X,td,V 


= (w + A x ^v 2 ) _°i^ + fi J dV 


(v - v') 2 G 


tjX^'jV' 


(2.16) 


where we have included a possible spatial dependence of the electron density in the form of 
A x depending on location in the 2D space. The neutrino velocity domain of integration is 
determined by the dimensionality of the chosen problem and if multi-angle effects are to be 
considered. 


2.5 Mass ordering 

In a two-flavor system, one important parameter for matter effects in general and for self- 
induced flavor conversion in particular is the mass ordering. In our context the question is 
if the dominant mass component of v e is the heavier one (inverted ordering) or the lighter 
one (normal ordering). Traditionally “mass ordering” is also termed “mass hierarchy” and 
we denote the two cases as IH (inverted hierarchy) and NH (normal hierarchy). We are 
concerned with 1-3-mixing, the corresponding mixing angle is not large, and so it is clear 
what we mean with the “dominant mass component.” 

Our equations are formulated such that they apply to IH, the traditional case where 
self-induced flavor conversion is important in the form of the bimodal instability. Of course, 
it has become clear that NH is actually the more interesting case. For NH, Am 2 is negative, 
but we prefer to consider Am 2 a positive parameter. Therefore, NH is achieved by including 
explicitly a minus sign on the r.h.s. of equation (2.10). This change of sign translates into a 
minus sign for u in the first bracket in equation (2.16). 

For flavor conversion, it is irrelevant if neutrinos oscillated “clockwise” or “counter 
clockwise” in flavor space, i.e., in equation (2.15) we may change i —> — i or H — > — H without 
changing physical results. However, the relative sign between oj and A and /i is crucial. 
Therefore, switching the hierarchy is achieved by 

IH —» NH: A —>■ —A and p-t-fi. (2.17) 

In our stability analysis we will consider the parameter range — oo < n < +oo and —oo < 
A < Too as these are simply formal mathematical parameters. Physically both parameters 
being positive corresponds to IH, whereas the quadrant of both parameters being negative 
corresponds to NH. 

2.6 Linearization 

As a next step, we linearize the EoMs in the sense that the complex off-diagonal element of 
every G is supposed to be very small compared to its diagonal part. We write these matrices 
explicitly as 
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(2.18) 




where g is a real and G a complex number and all quantities carry indices (i,x,w,v). To 
linear order in G we then find the EoMs 


i(d t + v • V x )5i iXja; ,v 

i(d t + v • V x ) G t , x , w ,v 



G 




(2.19a) 


(2.19b) 


Up to normalization, the “spectrum” gt,x,u,v is essentially the phase-space density of all 
neutrinos. It is not affected by flavor conversion, but evolves by free-streaming if it is not 
homogeneous. 


2.7 Homogeneous neutrino and electron densities 

We are primarily interested in self-induced instabilities. Disturbances in the neutrino density 
and/or the electron density will certainly exist in a real SN and can play the role of seeds 
for growing modes. However, if these disturbances are small, it is unlikely that they will 
be responsible for instabilities themselves. Henceforth we will assume that the neutrino and 
electron densities do not depend on the transverse coordinate x, although the flavor content 
may well depend on x. Free streaming does not change the density if it is uniform. As a 
consequence, gt,x,u,v does not depend on x or t and likewise, A x does not depend on x. 

With this assumption, g W)V describes the initially prepared neutrino distribution, i.e., 
their density in the phase space spanned by co and v. Inspecting the first integral in equa¬ 
tion (2.19b), we may write the three independent terms as 


yoj,v 


( 2 . 20 ) 


Here, e represents the “asymmetry” between neutrinos and antineutrinos. The second term, 
ei, represents a neutrino current which exists if their distribution is not isotropic and not 
symmetric between neutrinos and antineutrinos. Overall, the first term in square brackets 
becomes 

co + ^Av 2 + ^e/rv 2 — fie i • v + \e 2 g ■ (2-21) 

The last term is simply a constant and can be removed by changing the overall frequency of 
the rotating frame. Defining 

A = A + eg (2.22) 

the term in square brackets effectively becomes co + ^Av 2 -gei • v. We may also return to 
the notation used in our previous papers and define 


S-, 


G 


t,X.,LU,V — 


t,x,o;,v 


(2.23) 


The linearized EoM then takes on the more familiar form 


i(dt T v • V x ) St y x.,ui,v — (co T 2 Av ■ v) St,x,u,v j' dT 2^ ^ ) 5 , a/,v , 'S , t,x,a/,v' • 

(2.24) 
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This equation corresponds to equation (6) of reference [45]. Besides the streaming term 
(the gradient term on the l.h.s.) that we have now included to deal with self-induced in¬ 
homogeneities, we have also found the additional term g,e\ ■ v which is unavoidable in a 
non-isotropic system, irrespective of the question of homogeneity. This neutrino flux term is 
missing in reference [45]. The presence of this term modifies the eigenvalue equation for a 
non-isotropic system. 

2.8 Spatial Fourier transform 

We can now perform the spatial Fourier transform of our linearized EoM of equation (2.24). 
It simply amounts to replacing the spatial dependence on x of S by it dependence on the 
wave vector k and v • V x — > iv • k, leading to 

iSt, k,w,v = (w + pv 2 + k • v) S tikiWiV - n duj' dV |(v - v') 2 g^ySt^yy , (2.25) 


where k = k — ge\. Therefore, the wave vector k has the same effect as a neutrino current. 
Including an electron current would appear in a similar way. However, in the following studies 
of explicit cases we will focus on isotropic distributions and worry primarily about self-induced 
anisotropies, not about the modifications caused by initially prepared anisotropies. 

In equation (2.25) we have spelled out the meaning of the phase-space integral f dT' = 
f dui' f dv. Notice that the meaning of j dY' has changed in the course of changing variables 
that describe the neutrino modes. All phase-space factors and Jacobians have been absorbed 
in the definition of the effective neutrino-neutrino interaction strength g as well as the nor¬ 
malization of the “spectrum” g,^, v . In particular, if we begin with an ensemble consisting of 
only v e and v e and no v x or u x , then our normalizations mean that du f dv g UiV = 1. In 
this latter integral, we have only included positive frequencies (neutrinos, no antineutrinos) 
so that this normalization coincides with our definition that g is normalized to n Uf; . 

2.9 Oscillation eigenmodes 

In order to find unstable modes we seek solutions of our linearized EoM of the form k w v = 
leading to an EoM in frequency space of the form 

(pv 2 + k- V + cu - ft) = g J dJ J dv' |(v — v') 2 g^'y Qn,kyy ■ (2.26) 

Eigenvalues Q = 7 + in with a positive imaginary part represent unstable modes with the 
growth rate k. 


2.10 Monochromatic and isotropic neutrino distribution 

In our explicit examples we will always consider monochromatic neutrinos with some fixed 
energy, implying a spectrum of two oscillation frequencies oj = Two- Assuming that the 
energy and velocity distribution factorize, we may write the spectrum in the form 


fih>,V — h(jj f v . 

The monochromatic energy spectrum is 

huj = —a5{uj + cu 0 ) + 6(io - w 0 ), 


(2.27) 


(2.28) 


- 11 - 


meaning that we have a antineutrinos (frequency w = —wo) for every neutrino (w = wo). 
The spectral asymmetry is e = 1 — a. 

We will consider isotropic velocity distributions which, in addition, are uniform, cor¬ 
responding to blackbody-like angular emission in the SN context. In this case, / v = l/r v , 
where T v is the volume of the velocity phase space. The eigenvalue equation (2.26) finally 
simplifies to the form in which we will use it, 


(5 A V 2 + k • V -|- w - D) Qn, k,Lj, v = At J dJ h,j J dV (v - v ') 2 Qn,k,<y,v' • (2-29) 


We now consider systematically different cases of velocity distributions. 


3 One-dimensional system 

As a first case study we consider a ID system, i.e., the toy model of “colliding beams” that 
has been used in the recent literature as a simple case where one can easily see the impact of 
spontaneous spatial symmetry breaking [46, 51-53]. We go beyond previous studies in that 
we include the multi-angle matter effect and study the “footprint” of the various instabilities 
in the two-dimensional parameter space —00 < n < +00 and —00 < A < +00. This schematic 
study already leads to the conclusion that essentially the largest-scale instabilities are “most 
dangerous” in the context of SN neutrino flavor conversion. 

3.1 Single angle (v = ±1) 

3.1.1 Eigenvalue equation 

We begin with ID systems, i.e., colliding beams of neutrinos and antineutrinos with different 
velocity distributions. The first case is what we call “single angle,” a nomenclature which 
refers to the zenith-angle distribution of SN neutrinos. As we have explained, in our way 
of writing the equations, “single angle” means that the neutrino velocity distribution has 
|v| = 1. In our first ID case this means we consider two colliding beams with v = ±1. 
Matter effects can be rotated away. 

The eigenfunction Qn k,u,v now consists of four discrete components. We denote these 
four amplitudes with the complex numbers R for right-moving (v = +1) neutrinos (w = 
+wo), R for right-moving antineutrinos, and analogous L and L for left movers. Our master 
equation (2.29) then reads 

/wo + k 0 ~n not \ 

0 —wo + k —\i not 

—fj, fia wo — k 0 

\ — n fxa 0 —ujQ — k) 

corresponding to the equivalent result of Duan and Shalgar [52]. The eigenvalues D are found 
from equating the determinant of the matrix in square brackets with zero. This condition 
can be written in the form 

1 a \ ( 1 a 

—k + u 0 — D —k — wo — D ) \k + ujq — D k — wo — D 

This expression depends only on /.i 2 and therefore yields identical eigenvalues for positive and 
negative //, i.e., for both neutrino mass hierarchies, as noted by Duan and Shalgar. It is also 
even under k —> — k as it must because the system was set up isotropically, so the eigenvalues 
cannot depend on the orientation of k. 



[R\ 

R 

L 

W 


= 0 . 


(3.1) 
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3.1.2 Homogeneous mode (k = 0) 

For the homogeneous mode, k = 0, this eigenvalue equation simplifies considerably. We 
already know that we have the same solution for positive and negative //, where the latter is 
the left-right symmetry breaking solution discovered in reference [46] . We here limit ourselves 
to /t > 0 and need to solve the quadratic equation 

wg - H 2 - /i[(l + a)w 0 + (1 - a)H] = 0 . (3.3) 


It has the solutions 


(l-q)„ 

2 


wo + 


(1 — a) n 


— 2 /two . 


(3.4) 


Unstable solutions exist for 

2 J/_ 2 . . 

(1 + ,/S) 2 wo (1 - yfa/ 1 ’ 1 ; 

which for a = 1/2 is the range 12 —8\/2 < h/ojq < 12 + 8y/2 or numerically 0.6863 < /t/wo < 
23.31. The maximum growth rate is 


^max 


2 y/a 
1 — a 


w 0 • 


(3.6) 


For a = 1/2 this is K m ax/wo = 2\J~2 ~ 2.828. The maximum growth rate occurs at the 
interaction strength 

_ 2(1 +a) 

/bwc— (1 _ a ) 2 - v-n 

We show the growth rate normalized to its maximum in figure 2 as a function of /i//hc max . 
For this normalization, the unstable range is 


1 _ 2yjfa < n < 1 2\Jkx 
1 + a ^K max 1 + a 


(3.8) 


If we use a = 1 — e and expand to lowest order in e, this range is e 2 /8 < /j//i Kmaj( < 2 — e 2 / 8 . 
Therefore, even if e is not very small (e = 1/2 in figure 2), the unstable range is close to its 
maximum range from 0 to 2 . 


3.1.3 Inhomogeneous modes (k > 0) 

The quartic eigenvalue equation (3.2) is not easy to disentangle. However, for large k it sim¬ 
plifies and can be solved. We may guess that, for large k, the real part of H is approximately 
— k and, without loss of generality, we may go to a rotating frame such that H H — k. 
Moreover, based on numerical studies, Duan and Shalgar [52] have conjectured that for large 
k, the unstable //-range scales with y/k. This observation motivates us to write /i = my/cook 
and, without loss of generality, the eigenvalue equation reads 

wo awo \ / k ak 

o/q — H — wq — ClJ \2k + loq — H 2k — wq — H 
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Figure 2. Growth rate n for the unstable mode in the homogeneous (k = 0) ID case for a = 1/2. 
The maximum growth rate K max is given in equation (3.6), the corresponding interaction strength 
in equation (3.11). 


Now we can take the limit k —> oo and approximate the second bracket as (1 — a)/ 2. The 
resulting quadratic equation is now easily solved and yields 


Q e 2 m 2 ^ y/16 — 8m 2 (2 — e)e + e 4 m 4 
loq 4 4 


(3.10) 


One can try the same exercise with the opposite rotating frame = 17 + k and finds a 
similar-looking result, where however the argument of the square-root is always positive, i.e., 
there is no unstable mode with this property. 

From equation (3.2) one finds that the maximum growth rate n mSLX is the same as in 
the homogeneous case of equation (3.6), i.e., for large k the maximum growth rate does not 
depend on k and is the same as for k = 0. However, it now occurs at the interaction strength 


l? K 

r lx n 


4 (1 + a) 
(1 — a) 3 


cook. 


(3.11) 


The unstable range is the same as given in equation (3.8) if we substitute /V/bt max with 
(ju/^ Kmax ) 2 . The growth rate as a function of [i is the same as shown in figure 2 if we 
interpret the horizontal axis as (/i//i Kmax ) 2 with our new /r Kmax . 

For intermediate values of k, the maximum growth rate deviates slightly from the two 
extreme cases. It is somewhat surprising that the structure of the eigenvalue equation is such 
that the maximum possible growth rate depends only on the vacuum oscillation frequency 
too and a, but not on the potentially large scale k. 


3.2 Multi-angle effects (0 < v < 1) 

3.2.1 Eigenvalue equation 

We may now study the multi-angle impact of matter, in the spirit of the SN system, in our 
ID model by extending the velocity integration over the entire interval — 1 < v < 1 instead 
of using only the modes v = ±1. One approach is to introduce n v discrete velocities on 
the interval 0 < v < 1, i.e., 2 n v modes on the interval — 1 < v < 1, leading to a total of 
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4 n v discrete equations. One may then proceed to find numerically the eigenmodes. This 
approach is marred by the appearance of spurious instabilities and one may need a large 
number of modes to obtain physical results [43]. 

Therefore, we usually avoid discrete velocities and represent the eigenfunctions k,^,v 
as continuous functions of their variables. Equation (2.29) reads for our specific case 

/ +OO y /* + l 

duj h^f — J dv (v v ) QQ,k,uj',v' • (3.12) 


Since our system was prepared “isotropic” in the sense of left-right symmetry, the instabilities 
cannot depend on the sign of k so that it is enough to consider 0 < k < oo. The r.h.s. of 
equation (3.12), as a function of v, has the form Aq + A\v + A 2 v 2 . Because 1, v. and v 2 are 
linearly independent functions on the interval —1 < v < +1, the l.h.s. must be of that form 
as well and we may use the ansatz 


Qfl,k,u,v — | 


Aq + AiV + A 2 V 2 


\\v 2 -\- kv + ui — If 


for the eigenfunctions. Inserting this form on both sides yields 


(3.13) 


A 0 + A]V + A 2 v 2 




dv' ( v 12 


2vv' + v 2 ) 


Aq + A\v' + A 2 v' 2 
\ A v' 2 + k v' + oo' — 12 


(3.14) 


This equation really consists of three linearly independent equations for the parts proportional 
to different powers of v, so we get three equations linear in Aq, A\ and A 2 . This set of linear 
equations is compactly written as 


( h I 3 I 4 \ 

1 - - 2 h - 2/ 2 - 2/ 3 

\ Io h h ) 



(3.15) 


where 


I 

1 n — - 


r+oo 


r+1 


duj h u 


dv 


\\v 2 + kv + u> — Q ’ 


(3.16) 


where we have dropped the prime from the integration variables. 

Notice that I n is odd under k —> — k if n is odd, and it is even if n is even. Moreover, in 
the determinant of the matrix in square brackets in equation (3.15), every term involving I n 
with an odd power of n involves another factor I m with rn odd. Therefore, the determinant is 
even under k —> —k, in agreement with our earlier statement that without loss of generality 
we may assume k > 0. 


3.2.2 Homogeneous mode (k = 0) without matter effects (A = 0) 

As a first simple example we consider homogeneous solutions (k = 0) in the absence of 
matter effects (A = 0). The latter assumption requires an exact cancellation A = A + = 0 

between the matter effects caused by the background medium and by neutrinos themselves. 
We consider this case only for mathematical convenience without physical motivation. The 
velocity integrals vanish for odd powers of v. For even powers, and using the monochromatic 
frequency spectrum of equation (2.28), we find 


In = 


I 1 


a 


+ 


1 


2(1 T n ) ycuo T 11 — 11 J 2(1 T tl) 




(1 + a) ojq — (1 — a) 11 


uj; 


-n 2 


(3.17) 
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and the eigenvalue equation corresponds to 


det 



(\ o |\ 

( 0 3 0 I [(1 + a) cuo + (1 

V 1 o' \) 


a) 0 ] 


= 0 . 


(3.18) 


Before searching for solutions, we may diagonalize the 3x3 matrix. This leads to three 
independent equations of the form of equation (3.3), where we need to substitute 




30 


x 


5 + 3\/5 >0, 

-10 < 0 , 

5-3\/5 <0. 


(3.19) 


Therefore, we get three instabilities: One for /x > 0, the usual bimodal instability (IH), and 
two negative-/! solutions (NH). The maximum growth rate is the same in every case as the 
one that was found in Sec. 3.1.2 and was given in equation (3.6). The exact unstable /x-ranges 
have changed according to the /x scaling provided by equation (3.19). For our usual example 
with a = 1 / 2 , the instability ranges for both hierarchies in the colliding-beam example were 
0.68 < |/x /ojq | < 23.31. After ^-integration they become explicitly 

1.76 < n/uo < 59.74, (3.20a) 

—69.94 < h/uiq < —2.06, (3.20b) 

-409.44 < n/uo < -12.05. (3.20c) 

We conclude that integrating over the velocity interval — 1 < v < +1 modifies the unstable /x- 
ranges, breaks the symmetry between normal and inverted hierarchy, and introduces another 
normal-hierarchy instability. 

Qualitatively, these results are analogous to the three types of instability discovered 
in the study of axial-symmetry breaking in the SN context [45]. The one inverted-hierarchy 
solution appearing in all cases is the bimodal instability and corresponds to the original flavor 
pendulum [21]. The first normal-hierarchy instability is what was termed the multi-azimuth 
angle (MAA) instability, although in our ID system we have only two “azimuth angles,” i.e., 
the two beam directions. The third instability, appearing in normal hierarchy for a much 
larger /x-range, is what was termed the “multi zenith angle” (MZA) instability. It requires, in 
the SN terminology, a nontrivial range of zenith angles, corresponding here to a non-trivial 
u-range, i.e., anything beyond the trivial v = ±1 velocity distribution. 

Notice that our ID MAA instability corresponds to an eigenfunction which is anti¬ 
symmetric in v (it breaks the left-right symmetry) and corresponds to the middle entry —2/3 
in the matrix of equation (3.18) which decouples from the remaining 2x2 block. The latter 
yields left-right symmetric solutions (even under v —> —v). In this sense, it is the MZA 
instability which corresponds, for normal hierarchy, to the bimodal solution. It exists only 
for at least three velocity modes, and of course requires the presence of two vacuum oscillation 
frequencies, here always chosen asw = ±u;o, i.e., we need at least a total of six modes, making 
a simple visual interpretation more difficult. 

The growth rates of all three intabilities as functions of /x are shown in figure 3 as blue, 
green, and orange lines, all of them having the same maximum. In the top panel, we overlay 
the two instability curves for the original colliding-beanr example, where v = ±1. Therefore, 
this panel directly illustrates the effect of the velocity integration (“multi-angle effects” in 
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Figure 3. Growth rate n for the unstable modes in the homogeneous (k = 0) ID case with a = 1/2, 
where fj, < 0 corresponds to normal and > 0 to inverted hierarchy. (Both k and /i are in units 
of u>q.) The colored curves common to all panels are the three instabilities which obtain after velocity 
integration (—1 < v < +1) with the instability ranges of equation (3.20). The overlaid black instability 
curves are for discrete velocity bins, where the number of bins n v = 1, 2 , and 10 as indicated in the 
panels. The overlay in the top panel ( n v = 1) is the colliding-beam example with v = ±1. With 
increasing n v (top to bottom), discrete velocity bins approximate a uniform distribution. 


SN terminology) in that the velocity integration takes us from the two black curves to the 
three colored ones. 

We have also studied the equations using a set of discrete velocities, where the v = ±1 
case is the simplest example with n v = 1 bins. (We count the number of bins in the range 
0 < v < 1 , i.e., there are equally many bins for negative velocities, and the total number 
doubles for our two frequencies u = ±wo-) Adding the intermediate values v = ±1/2 takes 
us to n v = 2, shown in the second panel of figure 3. It reveals that the symmetry between 
the hierarchies (/i —> — /i symmetry) is broken as soon as the velocity range is non-trivial 
and that there are two normal-hierarchy solutions. Increasing n v eventually approximates 
a uniform v distribution. A fairly small number of velocity bins is enough to achieve good 
agreement. We will see shortly that, including non-zero k and/or A, changes the picture 
because spurious instabilities appear. 
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3.2.3 Inhomogeneous modes (fc > 0) without matter effects (A = 0) 

Non-vanishing matter effects (A ^ 0) and non-vanishing inhomogeneities (k ^ 0) modify the 
eigenvalue equation in similar ways: The range of effective oscillation frequencies given by 
\\v 2 + kv + co increases considerably if k/uo 1 and/or A/cuo 1- Roughly one would 
suspect that significant collective phenomena require a neutrino-neutrino coupling exceeding 
this range of frequencies, i.e., a /r range exceeding something like the rms spread of this 
range. In this sense, one would expect that the /r-range of unstable solutions would be 
shifted roughly linearly with A and/or k. 

We first test this picture with A = 0 and k > 0. The u-integrals in equation (3.16) can 
be performed analytically (Appendix A) and the w-integration amounts to summing over two 
terms with uj = Two- However, finding the eigenvalues of equation (3.15) requires numerical 
tools. We use Mathematica and show the result in figure 4 for 0 < k/oj o < 100 as indicated 
above the curves. In the left panels, we only show the MZA mode. For relatively small k , the 
function shifts left and deforms somewhat, whereas for larger k values, it shifts left nearly 
linearly with k. We have checked that this linear behavior obtains numerically even for very 
large k values—we have tested values up to 10'. In other words, for large k , the instability 
curves are very similar as a function of fi/k, although they narrow somewhat. In contrast to 
the earlier two-beam example, there does not seem to be quite a universal function for large 
k. Moreover, in the earlier case, the scaling was with ^/y/uiok, i.e., the nontrivial v range 
has qualitatively changed the results with regard to the fc-scaling. 




Figure 4. Growth rates k for different wave numbers k as indicated above the curves. The other 
parameters are A = 0 and a = 1/2. The left panels show only the MZA mode, in the right panels we 
show the MAA mode (/i < 0) and the bimodal mode (n > 0). The effect of non-zero k is essentially 
to shift the curves and for large fc, the instability curves are similar as a function of \i/k. 
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In the right panels of figure 4, we show analogous results for the MAA mode (/i < 0) 
and the bimodal mode (/i > 0), which show analogous behavior. 

One can study the same case by solving the eigenvalue equation in terms of velocity 
bins, in analogy to what one would do if one were to solve the EoMs numerically instead 
of performing only a linear stability analysis. For k = 25 luq we show the growth rates for 
all solutions in figure 5 for different choices for the number n v of velocity bins. (We recall 



Interaction Strength p 


Figure 5. Growth rates for all modes that appear when we consider n v positive and n v negative-^ 
bins, using k = 25 ujq, A = 0, and a = 1/2. For growing n v , more and more spurious modes appear, 
but their growth rates become smaller and the physical modes begin to stick out. 
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that in our convention, the total number of positive and negative v bins is 2 n v .) The large 
number of spurious modes is a conspicuous feature of these plots, although for sufficiently 
large n v , the physical modes stick out. 

With non-vanishing k and/or A, the functional form of the eigenfunctions in equa¬ 
tion (3.13) no longer factorizes as a function of v and one of u. It is conceivable that 
spurious modes can be avoided, or their impact reduced, if one were to find a better way of 
discretizing the neutrino modes than by simple bins in velocity and frequency [44], 

It is noteworthy that for all modes, spurious or physical, the growth rates are of order 
WO) i.e., they do not inherit a larger frequency scale from k or, in later cases, from A. Even 
for huge values of k and A, this conclusion does not change and agrees with our explicit result 
in the two-beam example. 

3.3 Including matter (A ^ 0) 

Including matter in our “multi-zenith-angle” case has the effect of introducing both A and 
k in the denominator of the integrals of equation (3.16). They can still be done analytically 
(Appendix A), but lead to transcendental functions. Of course, numerically one can find the 
eigenvalues without much problem. The parameter A, like k , has the effect of broadening 
the effective range of oscillation frequencies and of shifting the unstable collective modes to 
larger values of |/x|. We study the homogeneous and inhomogeneous cases separately. 


3.3.1 Homogeneous mode (k = 0) 

For the homogeneous mode (k = 0), the eigenvalue equation (3.15) simplifies considerably 
because in this case I± = I 3 = 0 and we are left with 


(h 

0 

Ml 

1-0 

-2 h 

0 

. Vo 

0 

h) \ 



We need to solve the two equations 


(I '2 — l ) 2 = I0I4 and. I2 


- 1 / 2 . 


(3.21) 


(3.22) 


As an overview, we show in figure 6 a contour plot of the growth rate ac in the two-dimensional 
parameter space of the interaction strength /a and the effective matter density A = A + e/j,. 

As explained earlier, the first quadrant (/a > 0 and A > 0) corresponds physically to 
inverted mass ordering (IH), whereas the third quadrant (/a < 0 and A < 0) corresponds 
to normal mass ordering (NH). The other quadrants would be relevant, for example, for a 
background medium of antimatter. Mathematically, /jl and A are simply parameters which 
we leave unconstrained by physical considerations. As usual, we use a = 1/2 and therefore 
e = 1 — a = 1/2 so that matter-free space (A = 0) corresponds to the line A = /x/ 2 . 

For A = 0 , the growth rates as a function of /a were shown in figure 3 in the form 
of the colored curves. The effect of increasing |A| is to shift the unstable regions to larger 
values of |/x|, creating the butterfly image seen in figure 6 . For /i>0we obtain the bimodal 
instability which exists even in a single-angle treatment. For /x<0we have two multi-angle 
instabilities as indicated. The solutions shown in the upper panel derive from the first block 
in the eigenvalue equation (3.22), providing the bimodal and MZA instabilities. The MAA 
solution shown in the lower panel derives from the second block in equation (3.22). 
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Figure 6. Growth rate n of the ID instabilities as a function of p and A, assuming a = 1/2. 
Upper panel: The first block in the eigenvalue equation (3.22) yields the bimodal instability for 
p > 0 and the multi-zenith-angle (MZA) instability for p < 0. Lower panel: The second block in 
equation (3.22) provides the multi-azimuth-angle (MAA) instability for p < 0. Notice that the first 
quadrant (p. A > 0) represents IH, the third quadrant (p, A < 0) represents NH. 


The contours for the different instabilities are quite different. In the regime of large 
A, one can expand the eigenvalue equations in powers of 1/A and identify analytically the 
behaviour of the footprint of the instability regions in the /x-A-plane (Appendix C). We show 
the footprint of the contours of figure 6 on a logarithmic scale in Appendix C in figure 15 
together with the asymptotic large-A expansion. 

3.3.2 Inhomogeneous mode (k > 0) 

We next determine numerically the same instability footprints for non-vanishing wave num¬ 
bers k. We already know the impact of nonzero k for A = 0, i.e., the instability is shifted 
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to larger fi values. We do not expect a big difference for A k relative to the k = 0 case. 
These expectations are borne out by our results shown in figure 7. Considering first the 
simpler fj, > 0 half of the plot, we show the instability footprints for k = 0, 10 2 , 10 3 and 10 4 
as indicated in the plot. We can now easily diagnose the impact of small-scale instabilities: 
Essentially they fill in the entire space between the k = 0 footprints between the quadrant 
with positive and negative A, whereas the region above the k = 0 footprint in the upper 
quadrant, and the space below in the lower quadrant remains stable. The only small caveat 
is that in the upper quadrant, for k ~ A, there are “noses” of the footprint sticking into the 
previously stable region. So there is a narrow sliver of parameters above the k = 0 footprint, 
the envelope of the noses, which becomes unstable due to small-scale instabilities. 



Figure 7. Footprint of the ID instabilities (k > 10~ 2 ) in the /r-A-plane for a = 1/2 and the indicated 
values of k. The homogeneous case (k = 0) is the footprint of the contour plot of figure 6, here on a 
logarithmic scale. The corresponding large-A asymptotic results are shown in Appendix C, figure 15. 
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In the left half of the plot (// < 0) the situation is somewhat more complicated because 
of the presence of two instabilities. For large k, the footprint actually connects asymptotically 
to the k = 0 instabilities in a crossed-over way which we illustrate only by the k = 10 1 case. 
In other words, for large k, the MAA and MZA instabilites strongly mix with each other. 


4 Two-dimensional system 


We now turn to a 2D system, corresponding to the SN example where neutrinos propagate 
within the “expanding transverse sheet” moving outward. The radial motion is parameterized 
by our “time” variable, whereas “space” is represented by two transverse directions. This 
example corresponds, with properly scaled variables p, and A, to the usual treatment of self- 
induced flavor conversion in SNe, except that now we can include small-scale instabilities 
with nonvanishing wavenumber k. 


4.1 Single angle (|v| = 1) 

4.1.1 Eigenvalue equation 

We begin with the “single zenith angle case,” meaning that the rescaled neutrino speed within 
the transverse sheet is |v| = 1 and the matter effect can be rotated away. Our velocity phase 
space is the unit circle, described by an angle variable p which we can measure relative to 
k. As the system is initially prepared axially symmetric, all vectors k are equivalent—the 
eigenvalues depend only on k = |k|. The eigenvalue equation (2.29) therefore reads 

/ -bOO y r+TT 

duj h^f —— I d<*p (1 CpCpt SpSpf) Qq, : k, lu', ip' 5 ( 4 * 1 ) 

- OO «/ —7T 

where c v = cos p and = sin <p. 

Proceeding as in our earlier cases, we notice that the r.h.s. of equation (4.1) has the 
form A\ + A c cos p + A s sin p, i.e., a superposition of three linearly independent functions on 
the interval —it < p < +7r. Therefore, the l.h.s. must be of that form as well and we may 
use the eigenfunction ansatz 

^ _ A\ + A c c<p + A s S(p 

' c cQ,k,bJ,(p — j ' 

k Cp + (jJ — \ l 

We may insert this form on both sides, leading to three linearly independent equations, 
corresponding to the coefficients of the three functions 1, cosy? and siny?. We may write the 
three equations in compact form 


(4.2) 


l h 

Ic 

0 \ 

(M 


1 - \-I c 

dec 

0 

Uc = 0 , 

(4.3) 

. V 0 

0 

-W 

Us) 



where 


la — P 



dui hu 


1 

2 b 



fa(<P) 

k c<p + u) — Vt, 


(4.4) 


Here, fi(p) = 1, f c (p) = cos p, f cc (p) = cos 2 p, and f ss (p ) = sin 2 p. We have used that such 
integrals vanish if they involve a single power of sin p because this is anti-symmetric on the 
integration interval, explaining the zeroes in the matrix in equation (4.3). We also note that 
I ss = I\ — I cc , so we need only three different integrals. 
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4.1.2 Homogeneous mode (k = 0) 

We first consider homogeneous solutions (k = 0) where the angle integrals in equation (4.4) 
can be performed explicitly. Using the monochromatic frequency spectrum of equation (2.28), 
the eigenvalue equation becomes 


det 


'1 0 

ujq — Q 2 — p I 0 — ^ 


[(1 + a) wo + (1 — a) fi] 


0 -i 


= 0. 


(4.5) 


These are three independent quadratic equations of the now-familiar form of equation (3.3). 
The first line corresponds to the usual bimodal solution, the second and third line to two 
degenerate multi-azimuth-angle (MAA) solutions which are unstable for negative p (normal 
hierarchy). For these modes, the instability range is a factor of 2 larger. 


4.1.3 Inhomogeneous modes (fc > 0) 

For the inhomogeneous modes (k > 0), the required integrals entering the eigenvalue equation 
are of the form 


I = » 
a k 


/ -Too 

' 

-oo 


did hujF a 


to — n 
k 


(*+7r 


where F a (w ) = — 
27T 


dip 


fa(<P) 

cos p + w 


With our monochromatic spectrum equation (2.28) we arrive at 




I 

a k f “ V k > V k 

We define the auxiliary function of a complex argument w 

s(w) = \/w — ly/w + 1, 

which, for complex numbers, is in general not equal to Vw 2 — 1. We then find 


1 w w 2 

= —{ v j Fc = 1 — r > F cc = —w -| r 
s[w) s(w) s{w) 


(4.6) 


(4.7) 


(4.8) 


F ss = Fi - F cc = w - s(w ). (4.9) 


These expressions allow us to write the eigenvalue equations explicitly, involving only poly¬ 
nomials and square-root expressions. 

We now have three non-degenerate solutions, in contrast to the original ID example 
of Duan and Shalgar [52], one for positive p (inverted hierarchy) and two for negative p. 
We show contour plots of the growth rate n as a function of p and k for our usual example 
a = 1/2 in figure 8. The two quasi-symmetric regions correspond to the two solutions which 
correspond to those of the ID case, i.e., the left-right symmetric and anti-symmetric cases 
also shown in reference [52] in a similar plot. In addition we see a third solution which 
is genuinely a result of the spatial 2D geometry with non-vanishing wave-vector k. The 
eigenfunctions in this case are proportional to sin p where p is the angle between k and the 
velocity v of a given mode. 

The eigenvalue equation for the decoupled ss-block of equation (4.3), the new genuine 
2D solution, is explicitly 


(1 + a) w 0 - (1 


a 


) D + ay/—k — wo — Qy/k — uiq — D 

— \J—k T ujq — Q\/ k T idQ — D = 


e 


(4.10) 
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Figure 8. Growth rates for the 2D case with |v| = 1 (“single zenith angle”), using a = 1/2. The two 
quasi-symmetric regions correspond to the two instabilities which already appear in the ID case [52], 
although here the unstable /i range shifts with fc 3 / 4 . The third instability is genuinely 2D, it has no 
counterpart in the “colliding beam” examples, and its unstable fi -range shifts linearly with k. 


which is a quartic equation. It has unstable solutions for /x < 0. In analogy to the discussion 
of the ID case, we can extract the large-A; limiting solution with the substitution Q = fl — k 
and /x = —k/( 1 — a) — ^/4mwo^(l + «)/(l — a) 3 , where m is a dimensionless variable. The 
detailed factors are not crucial and are the result of some tinkering. After the substitution 
one expands equation (4.10) for large k and keeps only the term proportional to the highest 
power of k. One finds unstable solutions for 0 < m < 1 and the growth rate k/ujq = 
4y/m(l — m) a/(l — a 2 ) with K max /u ;o = 2ct/(l — a 2 ) which is 4/3 for our example a = 1/2. 
The unstable /x-range scales linearly with k. Notice that the asymptotic solution obtains 
only for very large A;-values because the dominant term in the expansion of equation (4.10) 
is proportional to k 3 / 2 , the next one proportional to k, and becomes relatively unimportant 
only very slowly. 

The 2x2 block in equation (4.3) leads to a much more complicated equation, but still 
consists of polynomials and square-roots. We follow a similar approach and substitute D = 
Q — k and (jl = m u/!/ 4 A; 3 / 4 /\/T — a, e xpand th e eigenva lue equation i n powers of k and keep 
only the largest term, leading to y/2 \/—ujo — &V ujq — & = m 2 —ljJq — D — ay/ ujq — D). 

This quartic equation provides the large-fc solution, which is symmetric for both hierarchies, 
i.e., symmetric under m —> —m. Because the second-largest power of k is only a power 1/4 
smaller than the largest, the asymptotic solution requires extremely large k-v alues. Notice 
that in the ID beam example, the unstable range scaled with A; 1 / 2 , while here it is A; 3 / 4 . 

4.2 Multi-angle effects (0 < v < 1) 

4.2.1 Eigenvalue equation 

We finally turn to the “multi zenith angle” case of transverse velocities within the full 2D disk 
described by |v| < 1, meaning that multi-angle matter effects are now included. We describe 
the velocity phase space by the speed v = |v| and an angle variable tp which we measure 
relative to k as before. Noting that (1/T V ) f dv = (l/i r) f// dtp f/dvv, the eigenvalue 


- 25 - 












equation (2.29) becomes 
( 2 ~\~ k V C(p T UJ 12 ) Qfl,k,uj,v,ip 


/ +°° r+K dp' /•! i ' ^ 9 ^ 

dJ hu> / — / duV [|(u /2 + u 2 ) - W (c^cy + s^v)] Qn,k,u',v’,<p', (4.11) 

-OO J-TT 71 Jo 


where Cp = cos p and sp = sin ip. The r.h.s. as a function of v and p is Aq + A 2 U 2 + A c v cp + 
A s v Sp, so the eigenfunctions are of the form 

n _ A 0 + A 2 v 2 + A c v Cp + A s v Sp 

n,k,uj,v,ip I\ v 2 + kv C p -\-w —12 

As usual, we insert this form on both sides, leading to four linearly independent equations, 
corresponding to the coefficients of the four functions 1, v 2 , vcp and v Sp. We may write the 
equations in compact form 



( 4 1 

il 

T c 

° ^ 



( A o\ 

1 - 

i\ 


12 

0 



a 2 

-2 i c 2 

- 2 1% 

-21 3 CC 

0 



A c 


V 0 

0 

0 

"2 If) 



W 


(4.13) 


where we have used that f dp' vanishes for integrals involving an odd power of sin p 1 . The 
integrals are now 


r +00 


In = /' 


duj h u 


r+n dp 
2ir 


dv 


fa{p) 


2 \v 2 + kv Cp + cj — 12 


(4.14) 


Here, f\{p) = 1, f c {p) = cos p, f cc {p ) = cos 2 p, and f ss (p ) = sin 2 p. We also note that 

rss T 1 TCC 

1 3 ~ 1 3 i 3 • 

4.2.2 Homogeneous mode (k = 0) with matter (A > 0) 

In the homogeneous case (k = 0), the function cos p in the denominator disappears, all angle 
integrals can be performed analytically, and I c n = 0. Therefore, equation (4.13) becomes 



(h h 0 0 \ 



Mo\ 

1 - 

h h 0 0 



a 2 

O 

1 

O 

O 



A c 


O 

O 

O 

1 

1 — 1 
CO 



\aJ 


= 0, 


where the integral expressions after performing the dp integration are 


r»+oo 


In — M 


duo hnjKn 


where K„ = dv -r-=- 


/ 0 


1 \v 2 + u — 12 


The velocity integrals are explicitly 


Ki = i log ( 


A 


k 3 = i 
3 A 


A- 5 = x 


1 - 


1 + 2(00 - 12) 

2 (oj - 12) 


A 

2(w - 12) 

A 


log 1 + 


A 


2(cj - 12) 


+ 


2(w - 12) 

A 


log 1 + 


A 


2(w - 12) 


(4.15) 


(4.16) 


(4.17a) 

(4.17b) 

(4.17c) 
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The previous ss and cc blocks in equation (4.13) have now decoupled from the rest and are 
degenerate, leading to the MAA instability. The remaining 2x2 block provides the bimodal 
and MZA instability. In other words, we now need to solve 

(/ 3 - l) 2 = hh and J 3 = -l, (4.18) 

in analogy to reference [45] with a slightly different notation. This entire development is very 
similar to the ID case. 

In the limit A —> 0, K \, and K§ approach 1/2, 1/4 and 1/6 times l/(w — D), which 
one can also find by setting A = 0 before doing the u-integrations. The matrix can be 



Figure 9. Growth rate n of the 2D instabilities as a function of p and A, assuming a = 1/2. 
Upper panel: The first block in the eigenvalue equation (4.18) yields the bimodal instability for /i > 0 
(inverted hierarchy) and the multi-zenith-angle (MZA) instability for p < 0 (normal hierarchy). Lower 
panel: The second block in equation (4.18) provides the multi-azimuth-angle (MAA) instability for 
p < 0. This figure is analogous to the corresponding ID case shown in figure 6. 
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diagonalized and we find three independent quadratic equations of the form of equation (3.3) 
where we need to substitute /r —> —///4 and /j, —> fx (3 ± 2\/3)/12. Following the same steps 
as in Sec. 3.2.2, the instability ranges for our usual example a = 1/2 are found to be 

1.27 < /x/wo < 43.28, (4.19a) 

-93.25 < n/io 0 < -2.75, (4.19b) 

-602.81 < /z/wo < -17.75. (4.19c) 


These results are numerically similar to equation (3.20) for the corresponding lD-case. 

In analogy to the butterfly diagram of the ID case (figure 6 ) we show in figure 9 a 
contour plot of the growth rate n in the two-dimensional parameter space of the interaction 
strength /i and the effective matter density A = A + e/r. The result looks qualitatively similar 
to the ID case. Again, for /r > 0 (inverted mass ordering), we obtain the bimodal instability 
while for fi < 0 (normal ordering) we find the MZA and MAA instabilities. 


4.2.3 Inhomogeneous mode (k > 0) without matter (A = 0) 

Next we consider the relatively simple case of k > 0 without matter. We may write the 
integrals of equation (4.14) in the form 


In = t [ dujh u A“, where 77“ 

k J —OO 




v n fa (y>) 

VC^ + W 


and we have introduced w = (ui — Q)/k. We find explicitly 


i v/1 - w^J-w(l + w) 

K 1 = w --, 

\JW 

, 2 W 3 yf — wV! — W 2 (1 + 2 w 2 ) 

K 3 = “V 4--> 

3 o\/w 

x _ 8 w 5 y/—wVl — w 2 (3 + 4u; 2 + 8ui 4 ) 
15 = T5~ + 15V^ 

K!} = ~ ~ VJ 2 — \J—W 2 \J\ — W 2 , 

1 2u> 4 + a/— tc 2 \/l — w 2 (1 + 2 w 2 ) 

= 4 3 ’ 

Ag G = —re — w 2 — \J—w 2 \/\ — w 2 ^j , 




rc(3 — 2tc 2 ) 



(4.20) 

(4.21a) 

(4.21b) 

(4.21c) 

(4.21d) 

(4.21e) 

(4.21f) 

(4.21g) 


Notice that Af + Af = A 3 1 . 

With the help of these analytic integrals it is relatively easy to solve the eigenvalue 
equation numerically. We show a contour plot of the growth rate k in the fi-k -plane in 
figure 10. The 3 x3 block in equation (4.13) provides three different solutions, i.e., one 
for /r > 0 (the usual bimodal solution in inverted mass ordering) and two solutions for 
fj, < 0 (normal ordering). The lxl block provides a further solution for ^ < 0. Figure 10 
corresponds to figure 8 in the single-angle case. In comparison, we have one more instability, 
now a total of four, of which three are for fi < 0 (normal mass ordering). 
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Figure 10. Growth rates for the 2D case with 0 < |v| = 1 (“multi zenith angle”), using a = 1/2. 
This figure is analogous to the single-angle case shown in figure 8, but now we have three instabilities 
for /r < 0 (normal mass ordering) and the usual bimodal one for fj, > 0 (inverted ordering). For all 
four instabilities, the unstable fj, range scales linearly with k as discussed in the text. 


One may also extract the large-fc asymptotic behavior (Appendix D). For the 3x3 block 
in equation (4.13) one finds that the system is unstable for 


fl — &i 


(^k + m 


for 0 < m < m mav i, where i = 1,..., 3 , 


(4.22) 


i.e., the unstable /r region scales linearly with k, in contrast to the corresponsing ID case. 
The values of the coefficients a* and m maX j j are given in Appendix D. 

For the lxl block in equation (4.13) one finds an instability on a very narrow strip 
around fi = —6 k. However, the maximum growth rate decreases with A; -1 / 2 so that in the 
limit k —> oo this instability disappears and we are left with those arising from the 3x3 block. 


4.2.4 Inhomogeneous mode (k > 0) with matter (A > 0) 

As a grand finale, we now turn to the most general 2D case with matter (A > 0) and 
inhomogeneities {k > 0). We need to find the zeroes of the determinant in equation (4.14) 
and write the integrals in the form 


/“ = = [ du>h u , where K“ 
A J — OO 


r +n ckp 

J-tt 2vr 



V n fa(<P) 

v 2 /2 + qv Cp + w ’ 


(4.23) 


where q = k/X and w = — These integrals can be performed analytically; we provide 

our results in Appendix A. 

We next determine numerically the instability footprints for non-vanishing wave num¬ 
bers k and show the result in figure 11, which looks qualitatively similar to the corresponding 
ID case that was shown in figure 7. For the simpler /r > 0 half of the plot, we show the 
instability footprints for k = 0, 10 2 and 10 3 as indicated in the plot. These k > 0 footprints 
fill the space between the k = 0 footprint and the horizontal axis. In addition, in the upper 
panel, there are small “noses” of the k > 0 footprints which slightly extend in the space 
above the k = 0 footprint, but this is a very small effect. 
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Figure 11. Footprint of the 2D instabilities (k > 10 -2 ) in the /z-A-plane for a = 1/2 and the indicated 
values of k. The homogeneous case (k = 0) is the footprint of the contour plot of figure 9, here on a 
logarithmic scale. The corresponding large-A asymptotic results are shown in Appendix E, figure 16. 


For [i < 0 the situation is more complicated because there are three instabilities. As 
in the ID case, for large k the footprint connects asymptotically to the k = 0 instabilities 
in a crossed-over way, which we illustrate for k = 10 3 . The main novelty of the 2D case is 
the appearance of another instability, which merges with one of the others when A> fc. In 
other words, one of the instabilities somewhat splits into two unstable ranges. We also recall 
that for very large k, the maximum growth rate of one of them decreases and vanishes for 
k —> oo, in which case we are back to a total of three instabilities. In the unphysical third 
quadrant, we notice somewhat pronounced “noses” of the k > 0 footprints. 

In all cases the main message is the same as in the ID case: The small-scale instabilities 
fill the space between the k = 0 instability and the horizontal axis, whereas the space between 
the k = 0 instability and the vertical axis remains stable. 
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5 Conclusions 


Several recent papers [46, 51-53] have studied the phenomenon of spatial spontaneous sym¬ 
metry breaking in the “colliding beam” model of neutrinos interacting with each other re- 
fractively. One important finding was that for any neutrino density (or in our nomenclature 
for any value of the neutrino-neutrino interaction energy fi) there is some range of spatial 
wave vectors k where the system is unstable with regard to self-induced flavor conversion. 
As a consequence, it seemed that an interacting neutrino gas would never be stable for any 
conditions, with potentially far-reaching consequences for SN physics. 

We have studied similar models, but including the multi-angle matter effect. We concur 
with the previous results in that smaller-scale modes are unstable at larger values of [i for a 
given matter density. This means that on a “footprint plot” such as figure 11, modes with 
k > 0 fill the space between the traditional footprint for k = 0 and the horizontal axis, but 
not the space towards the vertical axis. If we show the instability footprint in a plot like 
figure 1, adapted to more physical SN parameters, the large-fc modes extend the instability 
region in the direction of larger neutrino density for fixed matter density. 

Therefore, if the instability footprint of the traditional homogeneous (or rather spher¬ 
ically symmetric) mode does not intersect with the SN density profile, the large-fc modes 
are safe as well. In this sense, the traditional large-scale mode remains the most sensitive 
stability probe. On the other hand, if the physical SN profile of density and neutrino fluxes 
intersects any instability region, instabilities on a large range of spatial scales will occur and 
one would not expect any simple outcome of the flavor conversion process. 

Our analysis is based on a linearized stability analysis of a model which we have de¬ 
veloped in section 2. We have formulated the problem is such a way that it includes, as 
special cases, the “colliding beam” examples of the previous literature, allows the inclusion 
of multi-angle effects by a simple modification of phase-space integration, and formulates the 
SN case as a 2D system evolving in time, i.e., the non-trivial dynamical evolution is in the 2D 
expanding sheet of neutrinos as a function of SN distance. All cases of the previous literature 
are covered in a single simple formulation. 

In particular, our homogeneous, multi-angle, 2D system corresponds to the usual sce¬ 
nario described in the previous literature [45, 82], where both the zenith and azimuthal 
multi-angle effects are included. Note that the usual single zenith angle description of the 
early days of the collective oscillation discussions [22] cover only a small part of the parameter 
space, for example the bimodal instabilities shown in our figure 8. In addition, for the first 
time we have have considered a scenario where inhomogeneous modes are included in the 
description of SN neutrino evolution. We have found that the homogeneous mode remains 
the dominant source of instability. 

Still, in many ways our investigation is only a mathematical case study that may or may 
not apply to a realistic SN. Our main simplification is that we assume the flavor content of the 
SN neutrino stream to be stationary and to evolve only as a function of distance. Moreover, we 
assume a uniform boundary condition at the neutrino sphere, i.e., global spherical symmetry 
of neutrino emission. In other words, our case study still contains substantial and nontrivial 
simplifying assumptions to reduce the complexity of the full problem. Future work will have 
to go beyond some of these simplifications to develop a more realistic understanding of what 
really happens to neutrino flavor in the dense SN environment. 
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A Analytic integrals 


When searching for the complex eigenvalues D we need various integrals that can be found 
easily with Wolfram’s Mathematica. There can be issues about the validity of the analytic 
expressions in the complex plane, so we here give the integrals explicitly. 

In the ID case, for a non-vanishing k and in the absence of matter effects (A = 0), we 
need integrals of the form 

r+i v n 

fn{w)= / dv——, (A.l) 

7-1 v + w 

where w is a complex number. We first define two auxiliary functions 

L(w) = log (^-0 , (A.2a) 

A{w) = 2 + w [*7r sign Im(u;) — 2 arctanh(iu)] . (A.2b) 

The required integrals are found to be 


f 0 (w) = L(w ), (A.3a) 

fi(w) = A(w), (A.3b) 

fiiw) = —2 w + w 2 L(w), (A. 3 c ) 

h(w) = § + w 2 A(w) , (A.3d) 

fi{w) = — | (w + 3w 3 ) + w 4 L(w) . (A.3e) 


The actual argument will be of the form w = (oj — £l)/k. 

For non-vanishing matter effects (A ^ 0) and non-vanishing k, we need integrals of the 

form 


r+l 


dv 


9n(p,w ) = 

J_i v z + pv T w 

where w is a complex number and p is real. Again we define two auxiliary functions 
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The required integrals are found to be 


go{p, w) = 2 B(p, w ), (A.6a) 

gi(p,w) = -pB(p,w) + K(p,w) , (A.6b) 

g 2 (p,w) = 2 + (p 2 — 2w) B(p, w) — p K(p, w) , (A.6c) 

gz(p, w) = —2p — p (p 2 — 3 w) B(p, w) + ( p 2 — w) K(p, w) , (A.6d) 


g±(p, w) = | (l + 3 p 2 — 3 w) + ( p 4 — Ap 2 w + 2 w 2 ) B(p , to) — (p 3 — 2 pw) K(p, w) . (A.6e) 

The actual arguments are going to be p = 2k /A and w = 2(cu — 0)/A. 

In the 2D case to solve equation (4.14) we need the integrals defined in equation (4.21), 
i.e., integrals of the form 


K" 

- 1 - x n 


r +n dip 

J- tt 2vr 



V n fa{y) 

v 2 /2 + qv c v + w ’ 


(A.7) 


where = 1, fc(<p) = cosy;, / cc (y>) = cos 2 and / ss (y>) = sin 2 These integrals can be 

found analytically with the help of Mathematica. We first define two auxiliary functions 


A q , w = arctan 


q- — w 
V—w 2 


+ arctan 


1 - 2q 2 + 2w \ 

\/4 q 2 - (1 + 2 w) 2 J 


B q , w = - yv - (1 + 2u>) 2 

Our desired integrals are then found to be 
K 1 — A <? 

1 v 1 — j 

Ag = [Aq jU , + 2 (g 2 — tc) , 

—2\J—w 2 + (6g 2 — 6u; + l) B QtW + 4 (3g 4 — 6</ 2 u; + 2rc 2 ) A 9)tl 
1 {^Bq^yj + 2 q Aq tW ^ S w 


Kl = 


1 ^ 
2 


A- = 


K\ = 


1 + 


2 q 

2 \J—w 2 — (6q 2 — 2u> + l) Bq )W — 4 q 2 (3 q 2 — 4 w) A 


— 1 — Aw — 


AT = 


4q 

2V-w 2 - (6 q 2 + 2w + l) B q>w - 4 q 2 (3 q 2 - 2w) A q ^ 


S„ 


8 q 2 


(A.8a) 
(A.8b) 

(A.9a) 
(A.9b) 

(A.9c) 
(A.9d) 

(A.9e) 

, (A.9f) 


where S w = —isign(Imu>). 


B Frequently encountered eigenvalue equations 


Based on our “monochromatic” neutrino spectrum equation (2.28) with vacuum oscillation 
frequencies Two and a the number of antineutrinos relative to neutrinos, we constantly 
encounter eigenvalue equations of the form 


F 


2/2 

ojq — 0 



= 1 — a , 


(B.l) 
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where p is an interaction energy. In the simplest case, F(x) = x, this is the traditional 
eigenvalue equation for the flavor pendulum if we notice that our p = /r( 1 — a)/2 , whereas // 
is the traditional interaction energy. We also encounter F(x) = yfx and F(x) = log(x). To 
study these equations, we transform them to dimensionless variables by the substitutions 

_ 14" ol IT ol . . 

H = rruu o- and iI = wojq -, (B.2) 

1 — a l — a 


leading to 

( 2 m \ ( 2 m \ 

F V(l- a)/(l + a)-w) ~ aF {-(l-a)/(l + a)-wJ = 1 ~ a ' ( ‘ 3) 

These substitutions allow us to easily take the limit of a “symmetric” neutrino distribution 
with e = 1 — a —> 0. 

For the simplest case of the linear function F(x) = x, equation (B.3) becomes a 
quadratic equation. It has the solution 


w = —m ± \ — (2 — m)m T 


1 — a 
it ol 


(B.4) 


It has a nonvanishing imaginary part in the range 

( 1 - VS ) 2 „ .(1 + v ^) 2 

—--< m < —--. 

1 H - O. 1 “b Ol 


(B.5) 


We denote with Ki(a,m) the function which is the positive imaginary part of w. In the 
present case, it simply is 


iFi(a,m) = y(2-m)m- . (B.6) 

As a function of m, it is a semi-circle with center at m = 1 and radius 2y / a/(l T ct), i.e., 

A ™ = ( B .7) 

it ot 

This function is shown for several values of a in the top panel of figure 12. In the limit of 
equal neutrino and antineutrino densities, i.e., for a —> 1, the imaginary part is K(l,m) = 
— m)m. This limiting result can be found from equation (B.3) directly by substituting 
e = 1 — q, expanding the equation in powers of e, and keeping only the lowest power, leading 
to the equation 2m T 2 mw T w 2 = 0. 

To compare with Sec. 3.1.2, notice that p = fj,( l — a)/2 = mwo(lTa)/(l-a). Using this 
relationship between // and m as well as the one between Q and w reproduces the previous 
results. In particular, for a->lwe find the physical growth rate k = yj (2 fi — ojo)u>o which 
grows without limit for n —> oo. It is interesting that in terms of the scaled variable m 
one obtains, as a limiting results for a —> 1, a semi-circle for K(l,m) as a function of m. 
So in principle one could study all of our problems in the a —> 1 limit, and yet obtain 
representative results for a < 1, i.e., one could essentially eliminate the annoying parameter 
a and still obtain meaningful results for the asymmetric case. 
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Figure 12. The functions K(a,m) as defined in the text for the cases F(x) = x, x , and log(;r) as 
indicated, for a = 1, 0.3, 0.1, 0.03, and 0.01, in each case from outside in. 


The next more complicated case is the square-root function, F(x) = y/x, for which 
equation (B.3) takes on the form 


2 m 


(1 — a)/(l + a) — w 


1/2 


— a 


2 m 


— (1 — a)/(l + a) — w 


1/2 


= 1 — a. 


(B.8) 


It can be transformed to a quartic equation, but the results are too cumbersome to deal with 
and not particularly informative. Again we can expand this equation in powers of 1 — a and 
find, for the limit a —> 1, the cubic equation 2 w 3 + m(l + 2 w) 2 = 0. The imaginary part of 
the solution is 


K 1/2 (l,m) = 


8m (3 — 2m) + (2m) 2 / 3 [8m (9 — 4m) + 3\/81 — 48 m — 27] 


2/3 


4\/3 (2m) 1 / 3 [8m (9 - 4m) + 3^/81 - 48 m - 27] 


1/3 


(B.9) 


This function is shown in the second panel of figure 12 and looks almost like a semi-circle, 
but is not quite one. It is is nonzero for 0 < m < 27/16 = 1.6875 and takes on its maximum 

value of ye ^/|(69 + ll-v/3^) = 0.880086 at m = ^(3 + \/33) = 0.819803. In figure 12, we 
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show i^i/ 2 (a,m) also for several other value of a. The curves always begin at m = 0, i.e., 
there is no lower threshold in this case. 

We finally turn to the logarithmic case, F(x) = log(x), natural logarithm always un¬ 
derstood, for which equation (B.3) becomes 


( 2 m \ , f 2 m \ 

° g \ (1 — a)/(l + a) — w J _a ° g \ — (1 — ol)/( 1 + a) — w) = 1 " " ' ( ’ 10) 

There is no general analytic solution, but again we can consider the a —> 1 expansion where 
we find 1 = w[ 1 + log(— w/2m)\. The solution is w{m) = \/W{—e/2m), where e is Euler’s 
number and W(z) is the Lambert IT-function, i.e., W(z) is the solution of z = W e w . In 
MATHEMATICA it is implemented as IT (z) = ProductLog [z]. The positive imaginary part 
of our w(m ), i.e., K\ og (l,m) = Im[— 1/IT(—e/2m)], is nonzero for 0 < m < e 2 /2 = 3.69453 
and is shown in the bottom panel of figure 12. Again, it looks deceivingly like a semi-circle, 
but is not exactly one. Its maximum occurs at rn = 1.76684 and = 0.724611. 

We usually consider a = 1/2 as our main example. In this case, the eigenvalue equation 
can be solved analytically with the solution 


e — 3m ± W3m.(3m — 4e) 

^=V2 = -£- • 

The imaginary part is nonzero for 0 < m < 4e/3, has its maximum at m 
on a maximum value of 2/3. 


(B.ll) 
2e/3, and takes 


C Asymptotic solutions for ID and A —oo 

We can derive analytic asymptotic solutions for the ID case with matter, i.e., the large-A 
continuation of the contour plot figure 6. We begin with the bimodal and MZA instability for 
A > 0 and consider the first block of the eigenvalue equation (3.22). It is of the form l+C'i / u+ 
C2/r 2 = 0, where the coefficients C\ and C 2 depend on a, coo, D and A. We have evaluated the 
integrals according to the explicit transcendental functions given in equation (A.6). Inspired 
by numerical solutions, we assume that both the real and imaginary parts of the solutions fl 
remain of order ojq and do not become large as A — > 00, an assumption that later bears out 
to be consistent with the solutions. Therefore, we may expand Ci and C 2 in powers of A -1 
and find that the dominant terms are C\ oc A -1 and C 2 oc A~ 3 / 2 . In terms of a dimensionless 
interaction strength jl of order unity we write 

M= T ^(6A) 1 / 2 (2h 0 1/4 A S '' 4 , (C.l) 

where the exact coefficient was chosen for later convenience. The lowest-order term in C2/i 2 
no longer depends on A, whereas the lowest-order term in C\fi oc A -1 / 4 and slowly becomes 
small as A — > 00. To lowest order in A -1 , the eigenvalue equation is found to be 



This is an example of the type of equations that we always encounter in this context and 
which are discussed in Appendix B. 
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Figure 13. Growth rate n of the bimodal and MAA instabilities, using a = 1/2. The interaction 
strength is scaled according to equation (C.l). The blue curves show the asymptotic behavior for 
A —> oo, the other curves are for the indicated A values. 


The asymptotic solution derives from the term quadratic in /r and thus remains un¬ 
changed under fj, —> — //, i.e., it applies to both hierarchies. We show the asymptotic solution 
as a blue curve in figure 13 on a linear and logarithmic scale. We also show the growth rates 
for A = 10 2 , 10 4 and 10 6 where the solution is not symmetric under fj, —> —/j, because the 
linear term in // kicks in. We have already noted that one needs very large A values to obtain 
the asymptotic solution because the second-largest term only scales with A~ 1//4 relative to 
the dominant term. The asymptotic behavior is achieved for much smaller A values if // > 0. 
The growth rate vanishes completely above a certain |/x| value, but obtains nonzero values 
otherwise, i.e., there is no lower ft threshold. However, for ft < 0.5, the growth rate is a steep 
power-law of ft and can be taken to be effectively zero. 

For our usual example a = 1/2 we find numerically that the maximum growth rate 
occurs for ft = 1.494. Therefore, we find that 

H = ±4.911 Wq / 4 A 3/4 (C.3) 

gives us the locus of the maximum growth rate in the ^u- A plane for the bimodal and MAA 
solutions. The maximum value of ft before the growth rate becomes zero is 1.7724. On the 
small-/! side, the growth rate drops below k < 1/100, our usual criterion, at ft = 0.3478. 
Therefore, the footprint of the instability is the region between the lines ^ = 1.143 A 3 / 4 and 
5.826 A 3 / 4 , where both /r and A are given in units of the vacuum oscillation frequency ujq. 
This footprint is shown in the first quadrant (upper right) of figure 15. The corresponding 
footprint in the second quadrant (upper left) is also shown. 

We next turn to the MZA solution which exists only in inverted hierarchy (// < 0) and 
we consider the second block in equation (3.22). If we use ji = — A/[2(l — a)] the leading terms 
cancel, leaving us with a leading term of order A~ 1//2 . To obtain the lowest-order equation, 
we introduce another dimensionless parameter ft and write 


-A 


7r 


M = 


2 (1 - a) M 2 (1 - a) 2 

where, of course, the detailed coefficients in the second term are chosen for later convenience. 
One then finds a quadratic equation with solutions 


x/T 3 


a- 


\/woA, 


(C.4) 


f l 1 + a 2 — 2 ft 2 (1 + a 2 ) 

wq 1 — a 2 


± i 


4a 


1 — cr 


V 7 A 2 (i - A 2 ) • 


(C.5) 
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Notice that these solutions require 0 < p < 1 and we have always assumed 0 < a < 1. The 
imaginary part, as a function of p 2 , has the familiar semi-circular shape. In figure 14 we 
show it as a function of p (blue curve) and we also show the full solution for A = 10 3 and 
10 2 . The asymptotic solution is quite good for relatively small A values. In contrast to the 
other solutions, on a logarithmic scale the unstable range becomes very narrow as A —> oo. 



Figure 14. Asymptotic growth rate n for the MZA instability, using a = 1/2. The interaction 
strength is scaled according to equation (C.4). The blue curve shows the asymptotic behavior for 
A —> oo according to equation (C.5), the other curves are for A = 10 2 and 10 3 , from outside in. 

The maximum growth rate obtains for p = \/y/2. Therefore, for our usual example 
a = 1/2 we find that for the MZA solution, 


/r = —A 


7 T\ 2>UJq\/2 


(C. 6 ) 


gives us the locus of the maximum growth rate in the /i-A plane. The growth rate becomes 
exactly zero for p < 0 and p > 1, so the footprint (see upper-left quadrant in figure 14) is 
delimited by the curves /i = — A and /i = — A — The width of the footprint scales 

with V% i.e., on a logarithmic scale it becomes very narrow for large A. 

For A < 0 , the above approach does not lead to unstable solutions. Numerically we 
observe that for A —> — oo, the real part of the solutions approaches Re(f2) —> A/2, i.e., a large 
negative number. Therefore, to be able to expand the equation, we express Vt = \/2 + w ojq 
and seek self-consistent solutions with the dimensionless eigenvalue w of order unity. After 
expansion for A —> — oo, the asymptotic eigenvalue equations are 


log(l — w) — alog (—1 
1 — a 



(C.7) 


where 


a = log 


-|A) 

e A uj 0 ) 


1 

A 


or 


a = log 


-4A) 


3 + p 2 
(3 - p)p ’ 


(C. 8 ) 


where e is Euler’s number and as always the logarithm is with base e. This equation is one 
example for the type discussed in Appendix B. 

For our usual example a = 1/2, we can solve this equation analytically with the explicit 
result 


2 — e a ± iWe a (8 — e a ) 
w a =i/2 = -^- 


(C.9) 


- 38 











It has a nonzero imaginary part for —oo < a < log(8) = 2.0794, although it becomes 
exponentially small for a -C — 1. The maximum imaginary part obtains for a = log(4) = 
1.3863 and the maximum is 2. Therefore, the maximum growth rate obtains for 


A = 


or 




where 


A = - log(4) - 2 + log-= log 

wo / 


A = - 


2A\ 


3 + /r 
(3 - A)A ’ 


-A 

2 e 2 wr 


(C.10) 


(C.ll) 


Therefore, we have altogether three solutions, corresponding to the three instabilities, with 



Figure 15. Footprint of the ID instabilities in the /i-A plane for k = 0 (homogeneous mode) and 
a = 1/2 as explained in the text. The colored regions derive from a numerical solution, where the 
blue footprints correspond to the 2x2 block in equation (3.22), the red solutions to the lxl block. 
The grey regions show the asymptotic solutions in the large-A limit derived in this appendix. 
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(C.12a) 


maximum growth rates on the locus in the fi -A plane given by 


H = —2A 


3A+ v / 3(A + 2)(3A-2) 


I 1 — +2A a , 

, 0 \ 3A + y/3(A + 2) (3A — 2) 
M “ +2A 2/A^T) 


- 2A 1 


6A, 


(C.12b) 

(C.12c) 


where the limiting behavior is understood for A —> oo. Because A —>• —oo, the first solution 
corresponds to positive + and thus to the bimodal solution, the second and third solutions 
are the MAA and MZA instabilities, respectively. 

To draw the footprints in the lower quadrants of figure 15, we notice that k = 0 for 
a > log(8) and on the other side k < 1/100 for a < log(4 — \/39999/50) = —9.90348. 
Therefore, the asymptotic footprints are limited by 


ai = log(8) and 02 = log (4 — \/39999/50) (C.13) 


from which the limiting curves are extracted by solving equation (C.8) for fi. Once more we 
note that two of the footprints are “wide” and nearly symmetric between // —> — fi, whereas 
the third instability has a very narrow footprint. 


D Asymptotic solutions for 2D with A = 0 and k —> 00 


We are looking for the large- A; solutions of the 2D case without matter (A = 0 ). We need 
to find the zeroes of the determinant of the matrix in equation (4.13). We first look at the 
3x3 block and calculate it according to the explicit integrals that we have found. Next we 
substitute the variables as uj = 1, a = 1/2, D = —k + x, and 

fj, = a (k + mVk), (D-l) 


where a is a coefficient to be determined and overall the substitution for /i is an educated 
guess. Except for the choice a = 1/2, everything is still completely general. The unknown 
frequency to be found is x. Its imaginary part is the growth rate which we are looking for. The 
parameter m is an effective interaction strength because it gives us n in this parameterised 
form. 

Next we expand the determinant as a power series for large k and find to lowest nontrivial 

order 


det(3x3 block) = 


2880 - 480 a - 424 o 2 - 11 a 3 


2880 


3 i 

320 L 


V2a 2 {32 + a) (2y/x -1-y/x + l) —= +0(l/k). (D.2) 

J \ rZ 


For term proportional to 1 j\fk to dominate we demand the first term to vanish, giving us 
three possible values for a from the requirement 2880 — 480 a — 424 a 2 — 11a 3 = 0. The 
explicit results are quite complicated expressions. Numerically one finds 


ai 

= -37.1825, 

(D.3) 

02 

= -3.42115, 

(D.4) 

03 

= +2.05821. 

(D.5) 
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In other words, we have three asymptotic solutions, where one is for positive // and two for 
negative /j, as expected. 

If we now imagine that a is one of these solutions, the first term in the determinant 
vanishes and in the second term we can substitute a 3 = (2880 — 480 a — 424 a 2 )/ll to remove 
the a 3 term. In anticipation of the result we further introduce the quantity 


162^/3 [120 - a (20 + 3a)l 

777 — 

max 11 [1080 -a (120 + 53a)] ’ 

(D.6) 

which for our three possible a values are numerically 


TTlmax,1 = 1.23675, 

(D.7) 

m maXi 2 = 4.49396, 

(D.8) 

m max ,3 = 2.77208. 

(D.9) 

Then we are left with the equivalent of the determinant equation 


V^m = im ma , x (2\/x - 1 - Vx + l) . 

(D.10) 

It has the explicit solutions 


5 10 m 2 8y / m 2 (m 2 — m 2 iax ) 

3 Hm-niax ^W-max 

(DU) 


The solution has an imaginary part for 0 < m < m max - Therefore, the large-A: footprint of 
the three instabilities is limited by the lines 


H = a,ik 


and 


fl — Q>i 


(k + m max 



(D.12) 


For //-values between these lines, the system is unstable. 

Finally we turn to the lxl block in equation (4.13). We proceed with the same substi¬ 
tutions except for 

H = a{k + b ), (D.13) 


where for the moment we leave open what b is supposed to mean. Expanding the lxl block 
determinant in powers of large k, we here find 

det(l x 1 block) = - 1 — (—9 + b + 3x) — 


+ i 


. 2y/2( 


2(x - 1) 3/2 - (x + 1) 3/2 


Again we can get rid of the first term, this time by setting a = 
instability is for negative //. The remaining equation is 


Tj +0(l/t 2 ). (D.14) 

—6, i.e., the footprint of this 


det(lxl block) = (9 — b — 3x) \ — iAy/2 2(x — l) 3 / 2 — (x + 1) 3//2 

k 1 


J k 3 / 2 


+ 0(l/k 2 ). (D.15) 


The leading term does not provide an imaginary solution. In other words, for very large k we 
do not have an instability. If we keep both the leading and next to leading term, we finally 
need to solve the equation 


(9 — b — 3x) Vk 


14^2 


2(s - 


l) 3 / 2 - (x + 1) 3/2 


(D.16) 


Solving this equation actually leads to an asymptotic solution where the growth rate exists 
for a range of b- values. However, the maximum growth rate decreases with 1 /Vk. Therefore, 
we have overall four instabilities, but for k —> oo the one from the single block disappears. 
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E Asymptotic solutions for 2D with k — 0 and A —> oo 

We can derive asymptotic solutions for the 2D case with matter, i.e., the large-A solutions 
of the eigenvalue equation (4.15), corresponding to the two equations (4.18). We begin with 
the 2x2 block and A —> +oo. As in the ID case, we assume that D remains of order u>o, an 
assumption which is confirmed by the results. We express the interaction strength in terms 
of a dimensionless parameter £i in the form 

(E.l) 


l l = 


A* 


l — a 


A. 



Figure 16. Footprint of the 2D instabilities in the /i-A plane for k = 0 (homogeneous mode) and 
a = 1/2 as explained in the text. The colored regions derive from a numerical solution, where the 
blue footprints correspond to the 2x2 block in equation (4.18), the red solutions to the lxl block. 
The grey regions show the asymptotic solutions in the large-A limit derived in this appendix. 
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(E.2) 


To lowest order in A 1 the eigenvalue equation is, using w = 

log(l - w) - alog(-l - w) ( A A 

-= a , where a = log - 

l-a \2u 0 J 


^2/^Cb 

2(A- 1 ) 2 

a 2 


This result is identical with equation (C.7), but with a different expression for a. To 
draw the asymptotic footprints we simply need to solve for p using the limiting a-values 
given in equation (C.13). The result is shown in figure 16 as grey shaded regions in the 
upper panels, to be compared with the blue regions which derive from a numerical solution 
of the full eigenvalue equations. 

Next we turn to the lxl block for the limit A —> +oo and express the interaction 
strength in the form 


M = “ 


A + pu; 0 log(A/2o;o) 
1 — a 


(E.3) 


With p = 0 the eigenvalue equation is identically fulfilled to lowest order in A -1 , i.e., to lowest 
order unstable solutions require fi = —A/(1 — a). This simple behavior indeed corresponds 
to the very “thin” footprint shown in red in the upper left panel of figure 16. Including 
p ^ 0 leads to an approximate eigenvalue equation which is not very simple and does not 
lead to simple asymptotic solutions. Expressing /i in terms of p as in equation (E.3) we can 
numerially find the growth rate k as a function of p as shown in figure 17. It is clear that 
the instability footprint in the logarithmic figure 16 will be very narrow. We also notice that 
the maximum growth rate decreases with increasing A. (For all of the other instabilities and 
for a = 1/2, the maximum growth rate K max = 2u>o.) 

For the next cases we turn to the limit A — > — oo. In this limit, we write 11 = \/2 + wujo 
in analogy to the ID case. In the A —> —oo limit, the eigenvalue is characterized by w values 
of order unity. We also write the interaction strength again in the form of equation (E.l). 
The limiting eigenvalue equation is the same as in equation (E.2), but now with 


a = log 


A 

2wq 


2(A-i) 2 

(p - 4) p 


or 


a = log 



p + 1 
A 


(E.4) 


where the first expression applies to the 2x2 block, the second to the lxl block of the 
eigenvalue matrix. As before, to draw the asymptotic footprints we solve for jl using the 


£ 

15 


cr 


$ 

o 

I_ 

o 



Figure 17. Growth rate k for the instability deriving from the lxl block in equation (4.18), using 
a = 1/2, and solving the full eigenvalue equation. The interaction strength is scaled according to 
equation (E.3). The curves are for the indicated values of A. 
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limiting a-values given in equation (C.13). The result is shown in figure 16 as grey shaded 
regions in the lower panels, to be compared with the blue and red regions which derive from 
a numerical solution of the full equations. 
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